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Testing of RANS turbulence models for stratified 
flows based on DNS data 

By S. K. Venayagamoorthy f, J. R. Kosefff, J. H. Ferzigerf AND L. H. Shihf. 


1. Motivation and objectives 

Stably stratified flows such as those in the atmosphere or in large water bodies such as 
the ocean, lakes, estuaries and reservoirs are prevalent in the natural environment. The 
presence of the buoyancy force due to the stratification may have a substantial effect on 
the flow development and mixing processes, and hence influence the distribution of scalar 
substances such as pollutants, suspended sediments, in the environment. 

Today, there exist numerous turbulence models for calculating turbulent mixing in the 
environment. These models range from the simple eddy viscosity models to the more 
detailed large eddy simulations (LES) and direct numerical simulations (DNS) (Axell 
& Liungman 2001). However, both DNS and LES can be too computationally expen- 
sive (and often too idealized) for most geophysical and engineering applications. This 
limitation has restricted modelers to RANS approaches commonly based on turbulent 
kinetic energy (TKE) closure schemes. The most widely used RANS models today are 
two equation models which solve two transport equations for the properties of the turbu- 
lence from which the eddy viscosity can be computed. The best known of these models 
is the k-e model which requires the solutions of the turbulent kinetic energy equation (a 
component of essentially all current multi-equation models) and dissipation of turbulent 
kinetic energy equation (Ferziger et al. 2003). 

In most geophysical flows, turbulence occurs at the smallest scales and one of the two 
most important additional physical phenomena to account for is stratification (the other 
being rotation). In this paper, the main objective is to investigate proposed changes 
to RANS turbulence models which include the effects of stratification more explicitly. 
These proposed changes were developed using a DNS database on stratified and sheared 
homogenous turbulence developed by Shih et al. (2000) and are described more fully in 
Ferziger et al. (2003). The data generated by Shih et al. (2000) (hereinafter referred to as 
SKFR) are used to study the parameters in the k-e model as a function of the turbulent 
Froude number, Fru- A modified version of the standard k-e model based on the local 
turbulent Froude number is proposed. The proposed model is applied to a stratified open 
channel flow, a test case that differs significantly from the flows from which the modified 
parameters were derived. The turbulence modeling and results are discussed in the next 
two sections followed by suggestions for future work. 
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2 . Turbulence modeling based on DNS data 

2.1. The Data 

The DNS data of SKFR are used to develop parameterizations of modeling coefficients 
typically found in RANS models (e.g. C' /t ) as functions of quantities that define the local 
state of the turbulence ( Fr k etc.). 

SKFR performed DNS of stratified homogeneous turbulent shear flows. The data have 
been extensively discussed by SKFR and hence the discussion will not be repeated here. 
However, it is important to realize that the data provide all of the properties of the 
turbulence up to second order statistics. The subset of the data used in this study consists 
of the highest initial Taylor microscale Reynolds number runs (Re \ 0 = 89). The value 
of Re \ o = 89 is relatively high for direct numerical simulations but not high enough 
to produce results that are independent of Reynolds number effects. It is nevertheless, 
high enough that the effects of the other parameters should be accurately represented. 
All the data used for this study were taken from the latter parts of the SKFR runs in 
order to ensure that the turbulence was fully developed. The physical time was non- 
dimensionalized as St where S is the shear rate defined as: 


S = 


dU 

dz 


(2.1) 


As discussed by Ferziger et al. (2003), the turbulence does not become fully developed 
until sometime later than St > 2. Most of the runs were continued until St = 12-14. 
In our present study, only the data for times between St = 8 and the end of the run 
were used. To render the plots less confusing, the data is averaged over this time period. 
In total, 37 runs were used to derive the results given below. Each run is characterized 
by the initial Reynolds number (which is the same for all the runs) and the gradient 
Richardson Ri g (defined in equation 2.2) which has a fixed value for each run. 


N 2 

^ 9 ~ ~g2 

where N is the buoyancy frequency defined as: 


N = 


f-gdp\ 1/2 

V P dz) 


( 2 . 2 ) 


(2.3) 


2.2. k-e model 

k-e is a commonly used two equation model, of which many variations have been sug- 
gested. Here, we base our proposed modifications on the standard version of the k-e 
model. In the presence of stratification, the turbulent kinetic energy equation can be 
written as: 

Dk 

-^ = P~e-B + D k (2.4) 

where P is the rate of production of TKE, which is given, for simple shear flows like the 
SKFR flows, by: 

P = -tivys (2.5) 

B is the buoyancy flux given by: 

B = — — p'w ' . 

Po 


( 2 . 6 ) 
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C M 

C el 

Cf 2 

O'k 0~e 

0.09 

1.44 

1.92 

1.0 1.3 


Table 1. Values of constants in the k-t model (Rodi, 1980) 


Hfc is the transport term (equal to zero for the SKFR flows) modeled using the gradient- 
diffusion hypothesis as: 

*-£(££) <«> 
and e, the rate of dissipation of the turbulent kinetic energy, is modeled in an analogous 
way to equation (2.4): 


De Pe e 


2 c Be 

- W3-— 


D f 


( 2 . 8 ) 


where D e is the transport term (equal to zero for the SKFR flows) modeled again using 
the gradient-diffusion hypothesis as: 


8 _ 

dz 


D f = — [ -± — 


v t de 
a* dz 


(2.9) 


The eddy viscosity is then given by: 


v t = C^ — 


and the eddy diffusivity is: 


«t = 


vt 

Pr* 


(2.10) 


( 2 . 11 ) 


where Pry is the turbulent Prandtl number. The constants that are commonly used are 
given in Table 1. 

The value for the buoyancy parameter C& is a matter of much discussion. Various 
values have been suggested, e.g. Rodi (1987) suggests a value of 0 < C e 3 < 0.29, Baum 
& Caponi (1992) suggest C e 3 = 1.14 and Burchard & Baumert (1995) have argued that 
Cf 3 should have negative values. It turns out from our studies that this is a very sensitive 
parameter for stratified flows. 

2.3. Proposed parametrizations 

Several parameters have been suggested to characterize the effects of the stratification, 
the most obvious of which is the gradient Richardson number defined in equation (2.2). 
However, Ferziger et al. (2003) argue that this is not the best choice as it represents the 
forcing rather than the properties of the turbulence. As pointed out by SKFR, the tur- 
bulent Froude number gives better correlations. It can be defined based on the quantities 
computed from the k-e model as: 

, Fr ‘ = ^. <2 ' 12) 
Further, it is also noted that there is a gradient Richardson number at which turbulence 
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energy neither grows nor decays called the stationary Richardson number Ri gs . Holt et 
al. (1992) showed that Ri gs is a function of the Reynolds number at least at smaller 
values of Re. If, in equations (2.4) and (2.8), we set Dk/ Dt = De/Dt = 0, and eliminate 
e, then for homogeneous flows, the following relation is obtained: 



C e2 - C e i 

c e2 - c e3 


(2.13) 


where Rif s is the flux Richardson number. Applying the gradient diffusion modeling 
concept to both B and P in equation (2.13) yields: 


Ri 



(2.14) 


where it is found that Pry « 1 in the stationary flow cases. Hence we assume that Ri f s 
and Ri gs are equivalent in the discussion that follows. 

We choose to insist that C e3 be zero at the stationary state in order to satisfy the 
constraint that the model predicts the existence of a stationary state under the conditions 
determined by SKFR. Rearranging equation (2.13) and using the generally accepted value 
of 0.25 for the stationary Richardson number at high turbulence Reynolds number, we 
get: 


C e2 = 


Ce 1 

1 - Rifs 


(2.15) 


The functional dependence of the stationary Richardson number Ri gs on the turbulence 
Reynolds number Ret defined in equation (2.16) was given by SKFR as shown in equation 
(2.17). 



(2.16) 


Rlga 1 + 103/Pe fc 


(2.17) 


Any of the parameters in the model e equation can be allowed to vary as a function 
of the turbulent Froude number but only one parameter (or combination of parameters) 
can be derived from the model dissipation equation and the DNS data, enforcing a 
limitation on the choices that can be made. In this study, the parameters that are chosen 
to depend on the stratification are C e3 and C M , while C e i and C e2 are independent of 
the stratification. We do however include the effect of the Reynolds number on C e2 as 
discussed earlier through equations (2.15) and (2.17). 

The buoyancy parameter C e3 can be calculated from DNS data using the model e 
equation as a function of Frk- To accomplish this, we compute (1 /e)de/dt from the SKFR 
data by fitting a least square straight line to log e as a function of time over the time 
range used in all of the data fitting. We then solve equation (2.8) for C e3 using equation 

(2.15) for C e 2 (as discussed by Ferziger et al, 2003). The resulting values are plotted in 
Figure 1. Unfortunately, there is no clear trend in the data to suggest any definitive fit. 
However, the data indicates that C e3 is of order unity in the strongly stratified region 
and possibly shows a slight increase with increasing Fr &. The fit we used for this study 
is also shown in Figure 1 and is given by equation (2.18). This correlation is based on 
the observation that the mixing efficency peaks at about Frk of 0.4 - 0.5 for the SKFR 
data and C e3 should be zero there. It should be noted that the data does not extend into 
the very strongly stratified regime (i.e. very low Frk values). It is a regime that has a 
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Figure 1. Buoyancy parameter C € 3 as a function of the turbulent Froude number, o: values 
calculated from the SKFR data, : fit given by (2.18). 



- 


0 

0 

0 


0 

0 0 

0 


0 

10 

1 0 
1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

0 

\ 

8 


\o °Q0Pb 

x O OOqO / 
\ OO ^ ^ 


0 


- 


mixture of internal waves and turbulence and care has to be exercised in applying any 
turbulence model in this regime. 

! 1.44 for Fr k < 0.35 

1.44 - 9 .6{Fr k - 0.35) for 0.35 < Fr k < 0.50 , , 

6.4(Fr fe - 0.5) for 0.5 < Fr k < 0.80 ^ 

1.92 for Fr k > 0.80 


We note that Figure 1 indicates that a better fit to the data could be given but we felt 
it important to repeat the condition that C e 3 be zero at Fr k = 0.5. 

The eddy viscosity parameter C p obtained from the data is plotted as function of 
Fr k in Figure 2. The value of C /t obtained from the SKFR data at large Fr k (weak 
stratification) is lower than the typical value of 0.09. However, we fit the data such that 
C t , reaches the asymptotic value of 0.09 at high Fr k values. Thus: 

f 0.125 Frl + 0.014Fr fc for Fr k < 0.35 

C y, = < 0.006(Fr fe - 0.35)/(0.02 + 0.1(Fr k - 0.35)) + 0.02 for 0.35 < Fr k < 0.60 
[ 0.08 tanh(FYfc) + 0.01 for Fr k > 0.60 

(2.19) 

The scalar transport can be modeled using the turbulent Prandtl number defined in 
equation (2.11). We plotted the turbulent Prandtl number as a function of the turbulent 
Froude number as shown in Figure 3. The curve fit shown in the figure is: 


f 1.4 for Fr k < 0.35 

\ 1.4 - 0.55(1 - exp(-7(Fr fe - 0.35)) for Fr k > 0.35 


(2.20) 


where we have chosen to keep the Prandtl number constant in the highly stratified 
regions. 




132 


Venayagamoorthy, Koseff, Ferziger & Shih 



0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 


Figure 2. The eddy viscosity coefficient as a function of the turbulent Froude number, o: 
values derived from the SKFR data, : correlation (2.19) 



Figure 3. The turbulent Prandtl number as a function of the turbulent Froude number, o: 
values derived from the SKFR data, : correlation (2.20) 


3. Results and discussion 

In order to test the proposed parameterization and compare it with other models, it 
is essential that we apply it to a flow that is significantly different from the homogenous 
flows from which the parameters were derived. We chose the stratified open channel flow 
as our test case. This is a flow for which direct numerical simulations and large eddy 
simulations have been performed by Garg et al. (2000) and Shih (2003). 

The tests of the model were performed using a 1-D water column model called GOTM 
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(General Ocean Turbulence Model) developed by Burcharcl et al. (1999). The various 
turbulence closure schemes in GOTM are based on a modular format that enables easy 
incorporation in 3-D ocean circulation models and also allows for refinements or exten- 
sions to the turbulence models. The source code is available to the public on the Internet 
website www.gotm.net. We incorporated our proposed model into the GOTM code. Since 
none of the built in stability functions in GOTM modify their parameters based on depth, 
we had to adapt the code to modify the parameters in k-e model discussed above to be 
depth dependent based on the local turbulent Froude number. 

The test case we use is somewhat artificial. It is a pressure-gradient driven open channel 
flow in which the density is held fixed at both the lower solid boundary and the upper 
free surface and in some ways similar to the experiments done by Komori et al. (1983). 
The Reynolds number Re T and the Richardson number Ri T based on the shear velocity 
for all the test runs were 682 and 31 respectively. Previous studies by Garg et al. (2000) 
have shown that LES produces results that are in good agreement with DNS results for 
open channel flows. Hence we use an LES run with identical conditions to those in our 
test case to assess the predictions of our proposed model and that of the standard k-e 
model with constant stability functions. 

The flow was first allowed to develop to a converged solution (without stratification 
effects) using the standard k-e model after which the stratification was imposed. The 
velocity and density profiles at the initial state (i.e. after spin up) are shown in Figure 
4. Also shown in Figure 4, are the turbulent kinetic energy and dissipation profiles to- 
gether with the turbulent Froude number and the turbulent Reynolds number profiles. 
Superimposed on these figures are the instantaneous LES profiles. 

A total of four runs were done using different combinations of the modifications (out- 
lined above) to the k-e model so as to determine the most suitable combination of param- 
eters that can match as closely as possible the LES results. The test runs discussed in this 
paper are outlined in Table 2. Run 1 is based the standard k-e model with C e 3 = 1.44, 
while run 2 uses the proposed parametrizations. In run 3, we vary only the eddy viscosity 
parameter C ^ and turbulent Prandtl number Pr t and keep both C e 2 and C e 3 constant 
as in run 1. Run 4 is an arbitrary case that was chosen to investigate the effects of just 
varying the turbulent Prandtl number using the correlation given by equation (3.1). 

The velocity, density, turbulent kinetic energy, dissipation of the turbulent kinetic en- 
ergy, turbulent Froude number and turbulent Reynolds number profiles at non-dimensional 
times tu*/h = 4.1 and tu*/h = 14.6 for the standard k-e model (run 1) and the full 
modified version (run 2) together with the LES profiles are shown in Figures 5 and 6 
respectively. 

The full modified version of the proposed model (denoted by run 2) underperforms 
the k-e model in terms of predicting the velocity profile in the channel especially at later 
times, as seen in Figure 6, even though it appears to capture the free surface velocity 
better. However, it does significantly better in predicting the evolution of the density 
profile, especially during the transient stages of the flow. 

In run 3, both C e 2 and C e 3 are held at constant values equal to the ones used for the 
k-e test run 1. It can be seen from the density profile in Figure 7 that the mixing is 
now predicted better for the highly stratified regions (i.e. low Froude numbers) in this 
run. However, this case underperforms the standard k-e case for most of the lower half 
of the channel depth. The discrepancy between the LES and the RANS results given by 
runs 1 and 3 for the lower half of the channel depth indicates that in general the k-e 
model allows too much mixing to occur from the bottom boundary. There is a strong 




Figure 4. Initial profiles (a) velocity; (b) density; (c) tke; (d) dissipation; (e) turbulent Froude 
number Fru', (f) turbulent Reynolds number Rek- — : LES results, : k-e results. 


Run 

c „ 

c e2 

C e3 

Pr t 

Remarks 

1 

0.09 

1.92 

1.44 

0.85 

Standard k-e model 

2 

eq.(2.19) 

eq. (2.15) 

eq. (2.18) 

eq. (2.20) 

Full modified version 

3 

eq.(2.19) 

1.92 

1.44 

eq. (2.20) 

Varying C M and Pr t only 

4 

0.09 

1.92 

1.44 

eq. (3.1) 

Varying Prt only 


Table 2. Test runs used for evaluating the predictions by the modified k-e model 













Figure 5. Profiles at tu*/h = 4.1: (a) velocity; (b) density; (c) tke; (d) dissipation; (e) Fr k ; 
(f) Rek • — : LES results; : run 1; ••• : run 2. 


possibility that the prescribed asymptotic value of 0.85 that the SKFR data suggests 
for the turbulent Prandtl number is low even though there are good engineering data 
that suggest this value is reasonable. However, other observations such as those due to 
Hogstrom (1996) suggest a value of Pr t = 1. In run 4, we use a different correlation for 
Pr t such that it asymptotes to 1.0 at high Fr k (see equation 3.1). The results as shown 
in Figure 7 indicate closer agreement with the LES results compared to runs 1 and 3 
respectively. . 

Pr t = 0.4 exp(— 2. 5Fr k ) + 1.0 for all Fr k (3-1) 


4. Conclusions and future work 

Modifications of the k-e model to account for stratification based on the SKFR data 
have been suggested and tested in this study. The test runs done on the stratified open 










Figure 6. Profiles at tu*/h = 14.6: (a) velocity; (b) density; (c) tke; (d) dissipation; (e) Fry; 
(f) Rek ■ — : LES results; : run 1 ; ••• : run 2 . 


channel flow highlight the importance of correctly modeling the turbulent Prandtl number 
Pr t as well as the buoyancy parameter C e 3. The results suggest that the turbulent Prandtl 
number should be close to unity for the neutrally stable flows. Further, it appears that 
the buoyancy parameter C e 3 has to be prescribed as a value of the order of C e 1 in order 
to correctly model the effects of the buoyancy force in the dissipation equation. The 
lack of a clear trend in C e 3 as a function of Fr k from the DNS data suggests that the 
strategy of setting C e 3 equal to a constant is probably the most reasonable approach. The 
results of this paper are useful only in the regimes of weak to moderate stratification. 
Clearly further test cases should be performed where the effects of the stratification are 
more pronounced and where Fr k is small over larger fraction of the flow depth. This will 
allow us to evaluate the effects of varying C M as function of Fry more precisely. Further 
work is also required to capture the boundary layer dynamics properly so that momentum 
balance can be achieved properly. We are also unable to comment on the transport models 
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Figure 7. Velocity profiles (left) and density profiles (right). — : LES results; : runl; ■■■ : 

run 2; - • - : run 3; - • - in bold : run 4 
at tUt/h = 4.1 (top) and tu*/h = 14.6 (bottom). 


(which may also depend on stratification) as the SKFR database is for homogeneous flows. 
It is important to obtain data for strongly stratified and inhomogeneous flows and apply 
them to the development of better models for stably stratified flows. 
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